library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✔ ggplot2 3.2.1 ✔ purrr 0.3.3
## ✔ tibble 2.1.3 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
calidad_aire <- readRDS("data/02_estructura.RDS")
calidad_aire_orig <- readRDS("data/raw/calidad_aire.RDS")
estaciones <- readRDS("data/raw/coordenadas_estaciones.RDS")
El conjunto de datos calidad del aire contiene 128357 observaciones y 19 variables. En las secciones siguientes analizaremos los datos haciendo las operaciones de transformación y limpieza necesarias.
Existen 24 estaciones distintas cuya situación se puede ver en el siguiente mapa
library(leaflet)
leaflet(data = estaciones) %>%
addTiles() %>%
addMarkers(lng = ~longitud, lat = ~latitud,
label = ~as.character(estacion)
)
Podemos estudiar la calidad de los datos en las distintas estaciones. Por ejemplo, puede que no todas las estaciones sean capaces de medir todas las magnitudes. En el siguiente gráfico representamos el porcentaje de mediciones ausentes (o no válidas) para cada magnitud y estación
porc_na <- calidad_aire %>%
group_by(estacion) %>%
summarise_at(vars(so2:nmhc), function(x) mean(is.na(x)))
porc_na %>%
pivot_longer(
cols = so2:nmhc,
names_to = "magnitud"
) %>%
ggplot(
aes(x = magnitud,
y = as.factor(estacion),
fill = value)
) +
geom_tile() +
theme_minimal() +
labs(
title = "Porcentaje de mediciones ausentes para cada magnitud y estación",
fill = "",
x = "",
y = "Estación",
caption = "MSMK: taller de calidad del aire"
) +
scale_fill_continuous(trans = "reverse", labels = scales::percent, breaks = c(1, 0.5, 0), limits = c(1, 0)) +
theme(legend.position = "top")
En el gráfico anterior se puede ver que la calidad de la información de las magnitudes no, no2 y nox es muy alta. Sin embargo, para el resto de magnitudes no hay homogeneidad entre las estaciones. Además, en este proyecto buscamos predecir la calidad del aire para la ciudad de Madrid. Por lo tanto, la información deberá estar agregada de forma que tengamos un único valor por magnitud y día.
La pregunta es ¿cómo agregamos la información? La aproximación más habitual sería obtener el valor medio de cada magnitud en cada día. Ésta es la agregación que seguiremos en este proyecto. No obstante, la agregación es una cuestión importante en el diseño del objetivo del proyecto. Por ejemplo, el valor máximo podría ser otra forma de agregación, de forma que se pudiesen activar protocolos de contaminación en el caso de que en algún punto de la ciudad se sobrepasen ciertos límites y que, utilizando la media, podríamos diluir los valores extremos.
Agrupamos los datos con la media en el siguiente código:
calidad_aire <- calidad_aire %>%
group_by(fecha, anyo, mes, dia) %>%
summarise_at(vars(so2:nmhc), mean, na.rm = TRUE) %>%
ungroup()
Nos podemos hacer la pregunta de si la calidad de los datos ha ido evolucionando con el tiempo. A continuación representamos aquellos días para los que las magnitudes no están informadas aún habiendo agregado todas las estaciones
calidad_aire %>%
pivot_longer(
cols = so2:nmhc,
names_to = "magnitud"
) %>%
mutate(value = is.na(value)) %>%
filter(value == TRUE) %>%
ggplot(
aes(x = magnitud,
y = fecha,
fill = as.character(as.numeric(value))
)
) +
geom_tile() +
scale_fill_manual(values = c("1" = "black", "0" = "white")) +
theme_minimal() +
labs(
title = "Porcentaje de mediciones ausentes en cada estación",
fill = "",
x = "",
y = "Estación",
caption = "MSMK: taller de calidad del aire"
) +
theme(legend.position = "none")
Teniendo en cuenta que una de las variables que consideramos importante en nuestro estudio es pm25, esta variable no tiene una buena calidad en el periodo de los datos anterior al 2005. Por lo tanto, eliminamos aquellas observaciones anteriores al 1 de enero del 2005.
calidad_aire <- calidad_aire %>%
filter(fecha >= as.Date("2005-01-01"))
El conjunto de datos tiene una componente temporal importante. Debemos asegurarnos si están recogidos todos los días del periodo temporal.
Para ello, comenzamos ordenando el conjunto de datos por la variable fecha
calidad_aire <- calidad_aire %>%
arrange(fecha)
Y comprobamos que la diferencia entre dos días consecutivos es siempre 1
table(calidad_aire$fecha - lag(calidad_aire$fecha), useNA = "always")
##
## 1 <NA>
## 5508 1
Efectivamente, no hay días faltantes en los datos.
En la siguienta tabla se muestran estadísticos relevantes para cada variable.
skimr::skim(calidad_aire)
| Name | calidad_aire |
| Number of rows | 5509 |
| Number of columns | 18 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| numeric | 17 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| fecha | 0 | 1 | 2005-01-01 | 2020-01-31 | 2012-07-17 | 5509 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| anyo | 0 | 1 | 2012.05 | 4.35 | 2005.00 | 2008.00 | 2012.00 | 2016.00 | 2020.00 | ▇▆▆▆▅ |
| mes | 0 | 1 | 6.49 | 3.46 | 1.00 | 3.00 | 7.00 | 10.00 | 12.00 | ▇▅▅▅▇ |
| dia | 0 | 1 | 15.73 | 8.80 | 1.00 | 8.00 | 16.00 | 23.00 | 31.00 | ▇▇▇▇▆ |
| so2 | 0 | 1 | 8.17 | 4.31 | 1.60 | 5.40 | 7.40 | 9.56 | 36.38 | ▇▃▁▁▁ |
| co | 0 | 1 | 0.40 | 0.18 | 0.14 | 0.28 | 0.34 | 0.46 | 1.73 | ▇▂▁▁▁ |
| no | 0 | 1 | 26.72 | 31.38 | 1.58 | 8.08 | 14.42 | 31.04 | 259.69 | ▇▁▁▁▁ |
| no2 | 0 | 1 | 42.75 | 18.34 | 7.92 | 29.08 | 39.77 | 53.73 | 133.92 | ▆▇▃▁▁ |
| nox | 0 | 1 | 83.71 | 64.32 | 11.50 | 41.88 | 62.21 | 101.42 | 532.00 | ▇▂▁▁▁ |
| pm25 | 10 | 1 | 11.98 | 6.17 | 1.50 | 7.50 | 10.67 | 15.00 | 58.67 | ▇▃▁▁▁ |
| pm10 | 0 | 1 | 23.34 | 13.88 | 2.58 | 13.83 | 20.50 | 29.17 | 216.67 | ▇▁▁▁▁ |
| o3 | 0 | 1 | 46.23 | 22.46 | 2.57 | 28.00 | 47.79 | 63.64 | 110.07 | ▅▆▇▅▁ |
| tol | 20 | 1 | 3.15 | 2.43 | 0.20 | 1.57 | 2.50 | 3.93 | 31.20 | ▇▁▁▁▁ |
| ben | 19 | 1 | 0.64 | 0.44 | 0.12 | 0.35 | 0.52 | 0.78 | 8.60 | ▇▁▁▁▁ |
| ebe | 19 | 1 | 0.76 | 0.51 | 0.10 | 0.38 | 0.72 | 0.95 | 5.70 | ▇▁▁▁▁ |
| tch | 0 | 1 | 1.42 | 0.12 | 1.06 | 1.33 | 1.39 | 1.47 | 2.42 | ▃▇▁▁▁ |
| ch4 | 0 | 1 | 1.25 | 0.11 | 0.83 | 1.18 | 1.24 | 1.30 | 2.16 | ▁▇▁▁▁ |
| nmhc | 1 | 1 | 0.17 | 0.09 | 0.01 | 0.11 | 0.16 | 0.21 | 0.66 | ▆▇▁▁▁ |
En un proyecto real, habría que inspeccionar el comportamiento variable a variable (en este caso es asumible analizar 18 variables). Aquí solamente analizaremos la variable pm25.
pm25La distribución de la variable pm25 es la siguiente:
calidad_aire %>%
ggplot(
aes(x = pm25)
) +
geom_histogram() +
labs(
title = "Distribución de la variable pm2.5",
x = "",
y = "",
caption = "MSMK"
) +
theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 10 rows containing non-finite values (stat_bin).
La medición de la magnitud a lo largo del tiempo se puede ver en el siguiente gráfico
p <- calidad_aire %>%
ggplot(
aes(x = fecha,
y = pm25)
) +
geom_line() +
scale_x_date(date_breaks = "1 year", date_labels = "%Y")
p
La información del gráfico anterior es difícil de asimilar. En este caso podría ayudarnos si tuviésemos un gráfico interactivo. Podemos convertir un gráfico de ggplot2 a interactivo con ayuda del paquete plotly
plotly::ggplotly(p)
En el gráfico anterior se aprecia vagamente un comportamiento cíclico de forma que los valores parecen incrementarse en los meses centrales del año. Podemos analizar este comportamiento en el siguiente gráfico
calidad_aire %>%
ggplot(aes(x = mes, y = pm25)) +
geom_jitter(alpha = 0.25) +
geom_boxplot(aes(group = mes), alpha = 0, color = "firebrick", size = .75) +
scale_x_continuous(breaks = 1:12)
## Warning: Removed 10 rows containing non-finite values (stat_boxplot).
## Warning: Removed 10 rows containing missing values (geom_point).
A priori, el mes parece relevante en el comportamiento de la variable. En la fase de modelización comprobaremos si influye en la predicción o no.
Antes de modelizar es necesario abordar dos tareas:
El conocimiento que se tenga del contexto del proyecto en muchos casos es fundamental para obtener mejores modelos. Distinguimos dos tipos de generación de variables. El primero de ellos es la generación variables a partir de la información que contiene el propio conjunto de datos. El segundo será introducir información externa.
dia_sem y finde)Al tratarse de datos de calidad del aire, puede ser un factor relevante el día de la semana. A priori (y será algo que tengamos que comprobar en la fase de modelización), el mayor o menor uso de los vehículos, podría influir. Generamos la nueva variable dia_sem.
calidad_aire$dia_sem <-
lubridate::wday(calidad_aire$fecha,
label = TRUE,
week_start = 1)
Igual que hicimos anteriormente con el mes, podemos representar el comportamiento de la variables pm25 en función del día de la semana.
calidad_aire %>%
ggplot(aes(x = dia_sem, y = pm25)) +
geom_jitter(alpha = 0.25) +
geom_boxplot(aes(group = dia_sem), alpha = 0, color = "firebrick", size = .75)
## Warning: Removed 10 rows containing non-finite values (stat_boxplot).
## Warning: Removed 10 rows containing missing values (geom_point).
En el gráfico anterior, el sábado y el domingo están, en media, ligeramente por debajo de los otros días. Por este motivo, vamos a crear también una variable que indique si la observación se corresponde con el fin de semana.
calidad_aire$finde <- as.numeric(calidad_aire$dia_sem %in% c("sáb", "dom"))
*_med_mes)Aunque disponemos de la variable mes y el modelo de predicción la puede tener en cuenta, a veces ayuda a tener el valor medio de cada mes
calidad_aire <- calidad_aire %>%
group_by(mes) %>%
mutate_at(vars(so2:nmhc), list(med_mes = mean), na.rm = TRUE) %>%
ungroup()
*_inc1)Es habitual en el caso de series temporales incluir el incremento porcentual respecto del día anterior.
calidad_aire <- calidad_aire %>%
mutate_at(vars(so2:nmhc),
list(inc1 = function(x) (x - lag(x))/lag(x))
)
Se podrían generar muchas más variables y, conforme se va ganando experiencia en un determinado contexto, mejor sabremos qué variables suelen ser relevantes.
Recordemos que el objetivo es predecir el valor de pm25 en el día posterior. Eso significa que, si queremos predecir el valor para el 9 de febrero, tenemos que suponer que solamente conocemos la información hasta el día 8 de febrero. Por lo tanto, el valor que podremos utilizar de las variables para predecir un día, debe ser la información disponible hasta el día anterior. Por ejemplo, no podemos utilizar el valor de so2 en el mismo día que el valor de pm25 que queremos predecir. Necesitamos que el valor de cada variable esté retrasado en un día. Esto lo podemos hacer mediante la función lag. Para entenderlo, vamos hacer primero un ejemplo. Si tuviésemos el vector c(1,2,3,4), si aplicamos la función lag obtendríamos
lag(c(1,2,3,4))
## [1] NA 1 2 3
Obviamente, el valor anterior del primer elemento del vector es desconocido y por eso aparece como NA.
Ahora generamos una nueva variable _lag por cada variable que tengamos que retrasar. Por ejemplo
calidad_aire$so2_lag <- lag(calidad_aire$so2)
Este procedimiento habría que repetirlo demasiadas veces y sería demasiado pesado para hacerlo de forma manual. Para automatizarlo podemos utilizar la función mutate_at:
calidad_aire <- calidad_aire %>%
mutate_at(vars(so2:nmhc, so2_med_mes:nmhc_inc1), list(lag = lag))
Nota: lo que acabamos de hacer se puede traducir como: aplica la función lag a aquellas variables que están entre so2 y nmhc. Al utilizar list(lag = lag), cada variable que se crea termina en _lag.
Para finalizar la creación de estas variables, debemos eliminar todas las variables que no son lag y quedarnos solamente con pm25 que es la variable que queremos predecir.
calidad_aire <- calidad_aire %>%
select(fecha, anyo:dia, dia_sem, finde, ends_with("_lag"), pm25)
Es interesante conocer la correlación de la variable objetivo con respecto a las predictoras:
cor_pm25 <- cor(calidad_aire$pm25,
select(calidad_aire, ends_with("lag")),
use = "complete.obs"
)
sort(cor_pm25[1,])
## o3_lag o3_inc1_lag o3_med_mes_lag so2_med_mes_lag
## -0.26097887 -0.15548085 -0.09729124 0.09123486
## co_med_mes_lag pm10_med_mes_lag no2_med_mes_lag nmhc_med_mes_lag
## 0.10186458 0.10445794 0.11884661 0.11903891
## ebe_inc1_lag nox_med_mes_lag no_med_mes_lag ben_med_mes_lag
## 0.12204174 0.12592394 0.12718846 0.12758728
## ch4_med_mes_lag tch_med_mes_lag tol_med_mes_lag ebe_med_mes_lag
## 0.12866286 0.12872078 0.12938145 0.13110331
## ch4_inc1_lag nmhc_inc1_lag no_inc1_lag tol_inc1_lag
## 0.14396524 0.15538853 0.16101512 0.16662015
## no2_inc1_lag ch4_lag tch_inc1_lag so2_inc1_lag
## 0.17772042 0.17826189 0.18764322 0.19016680
## nox_inc1_lag pm25_med_mes_lag nmhc_lag ben_inc1_lag
## 0.19828819 0.20517579 0.20549780 0.22828950
## co_inc1_lag pm25_inc1_lag tch_lag pm10_inc1_lag
## 0.25573605 0.28546376 0.29525896 0.30514457
## ebe_lag so2_lag ben_lag tol_lag
## 0.36193726 0.39123372 0.39895123 0.48218583
## no_lag co_lag nox_lag no2_lag
## 0.49360402 0.51821069 0.52853479 0.55917953
## pm10_lag pm25_lag
## 0.61601447 0.69459147
Nota: en la correlación utilizamos
use = "complete.obs"para que no tenga en cuenta losNAen el cálculo.
Puedes ver que la mayor correlación de pm25 se da con pm25_lag. Podemos corroborarlo representando la relación entre ambas variables mediante un diagrama de dispersión:
calidad_aire %>%
ggplot(
aes(x = pm25_lag,
y = pm25
)
) +
geom_point(alpha = 0.25) +
geom_smooth() +
theme_minimal() +
labs(
title = "Relación de PM2.5 con el día anterior",
x = "PM25 del día anterior",
y = "",
caption = "MSMK: taller de calidad del aire"
)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 15 rows containing non-finite values (stat_smooth).
## Warning: Removed 15 rows containing missing values (geom_point).
Por último, como hacemos habitualmente, vamos a dividir el conjunto de datos en train y test. Entrenaremos con datos hasta 2019-09-01 y los restantes para test:
train <- calidad_aire[calidad_aire$fecha < as.Date("2019-09-01"),]
test <- calidad_aire[calidad_aire$fecha >= as.Date("2019-09-01"),]
Igual que hicimos en la fase anterior, almacenamos estos datos en el disco duro.
saveRDS(train, file = "data/train.RDS")
saveRDS(test, file = "data/test.RDS")